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ABSTRACT 

Many efforts have been made to model the mass distribution and dynamical evo- 
lution of the circumnuclear gas in active galactic nuclei (AGNs). However, chemical 
evolution is not included in detail in three-dimensional (3-D) hydrodynamic simula- 
tions. The X-ray radiation from the AGN can drive the gas chemistry and affect the 



-2- 



thermodynamics, as well as the excitation of the interstellar medium (ISM). Therefore, 
we estimate the effects (on chemical abundances and excitation) of X-ray irradiation 
by the AGN, for atomic and molecular gas in a 3-D hydrodynamic model of an AGN 
torus. We obtain the abundances of various species from an X-ray chemical model. A 
3-D radiative transfer code estimates the level populations, which result in line inten- 
sity maps. Predictions for the CO 7=1— >0to7 = 9— >8 lines indicate that mid-7 
CO lines are excellent probes of density and dynamics in the central (< 60 pc) region 
of the AGN, in contrast to the low-7 CO lines. Analysis of the Xco/o; conversion fac- 
tors shows that only the higher-7 CO lines can be used for gas mass determination in 
AGN tori. The [C ii] 158 //m emission traces mostly the hot (T/, > 1000 K) central 
region of the AGN torus. The [C ii] 158 /zm line will be useful for ALMA observations 
of high redshift (z > 1) AGNs. The spatial scales (> 0.25 pc) probed with our simula- 
tions match the size of the structures that ALMA will resolve in nearby (< 45 Mpc at 
0.01") galaxies. 

Subject headings: methods: numerical, galaxies: evolution, galaxies: ISM 



1. Introduction 



The formation and growth of a central black hole and its interaction with intense star-forming 
regions is one of the topics most debated in the context of galaxy evolution. There is observa- 
tional evidence for a common physical process from which most act i ve galactic nuclei (A GNs) 
and starbursts originate (e.g.. ISoltanI Il982l : iMagorrian & et al.l 1 19981 : iFerrarese & Merritti 12000, : 



Graham et al.ll200ll : iHaring & Rixll2004r) . A plausible scenario considers that starbursts, super- 
massive black hole growth, and the formation of red elliptical and submillimeter galax ies, are con- 
nected through an evolutionary sequence caused by rnergers betwe en gas-rich galaxies (|Hopkins et al. 
20061 12OO8I : iTacconi & etaDl2008l : iNaravanan et all 120091 UOm . Jn this scenario, the starbursts 
and (X-ray producing) AGN s seem to be co-eval , and the interaction processes between them 
(phase d and e in Figure 1 of iHopkins et al.ll2008|) . that dominate the formation and emission of 
molecular gas, is one of the long-standing issues concerning active galaxies. 

Numerous molecules tracing different (AGN and starburst driven) gas chemistry have been 
detected in Galactic and (active) extragalactic environments. Studies have shown that chem- 
ical differentiation observed within Galactic mo l ecular clouds is also seen at larger (~100 pc) 



scales in nearby galaxies (e.g.jHenkel et al.ll 19871: iNguyen-O-Rieu et al.lll99ll:lMartm et al.l2003 



Usero et al.l2004:lTacconi & et al.ll2008l : IPerez-Beaupuits et al 
van der Werf & et alboiO ) 



200ll2009Ll2010l : lBaan et al.l2010l : 
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The evolution of the ISM in the inner 100 pc regio n around a 10^ Mq) supermassive black hole 
(SMBH) was investigated by IWada & NormanI (|2002|) (hereafter |WN02|) using three-dimensional 
(3-D) Euler-grid hydrodynamic simulations. They took into account self-gravity of the gas, radia- 
tive cooling and heating due to supemovae (SNe). A clumpy and filamentary torus-like structure 
was found to be reproduced on a scale of tens of pc around the SMBH, with highly inhomogeneous 
ambient density and temperature, and turbulent velocity field. Their results indicated that AGNs 
could be obscured by the circumnuclear material. This represents theoretical supp ort for observa 



tional evidence showing that some AGNs are obscured by nuclear starbursts (e.g. jLevenson et al. 
200 li I2OO7I : iBallantvnelboOSi and references therein). 



Several efforts have been made to estimate the molecular line emission from the nuclear region 
in thes e 3-D hydrodynamic simu lations, an d to corn pare the results with observational data. For in- 



stance, [\AMa_&Tomis^a|(|2Q05[) (hereafter I WT05|) derived molecular line intensities emitted from 
the nuclear starburst region around a SMBH in an AGN. They used the 3-D hydrodyna mic sirn u- 
lations (density, temperature, and velocity field data) of the multi-phase gas modeled by |WN02| as 
input for 3-D non-LTE radiative transfer calculations of ^^CO and '^CO lines. They found that the 
CO-to-H2 conversion factor (X-factor) is not uniformly distributed in the central 100 pc and the X- 
factor for '^CO 7 = 1 — > is not constant with density, in contrast with the '^CO 7 = 3^2 line. 
Similar ly, the r ole of the HC N and HCO^ high-de nsity trace rs in the inhomogeneous molecular 
torus of lWN02l was studied by lYamada et al.l (120070 (hereafter I YWT07|) . These non-LTE radiative 
transfer calculations suggested a complicated excitation state of the rotational lines of HCN (with 
maser action) and HCO^, regardless of the spatially uniform chemical abundance assumed. 

However, all these previous efforts to estimate the molecular line emissions from the central 
100 pc of an AG N leave room for improvements. First of all, the radiative cooling in the simula- 
tions by lWN02l are not consistent with the chemical abundances in the cold and dense gas because 
(coUisional) formation and (radiative) destruction of H2 by far ultravi olet rad iation (FUV) was not 
included. Therefore, the cold and dense gas in the simulations by |WN02| does not necessarily 
represent the dusty molecular gas phase around an AGN. 

Hence, in order to study the distribut ion and structures of the v arious density regimes of the H? 
gas, th e 3-D hydr odynami c simulations of lWN02l were extended by lWada. Papadopoulos & Spaans 
(120091) (hereafter |WPSQ9|) to solve the nonequilibrium chemistry of hydrogen molecules along with 
the hydrodynamics. The formation of H2 on dust and its radiative destruction by far ultraviolet ra- 
diation (FUV) from massive stars are also included in the model by IWPS09l . This allows to track 
the evolution of molecular hydrogen and its interplay with the H i phase in the central 64x64x32 pc 
region. Thus, the radiative cooling in the model by WPS09I is more consis tent wit h the chemical 
abundances expected in the cold ISM, in comparison with the models by IWN02l Different SN 
rates and strengths of the uniform FUV field were also explored in order to study their effects on 
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the structures of molecular gas. 



On the other hand, the inhomogeneous density and temperature structures observed in the 3-D 
hydrodynamic models are not the only factors that drive molecular abundances and excitation con- 
ditions of molecular lines. There is observational and theoretical evidence in the literature that sup- 
ports different chemical evolution scenarios due to X-ray and UV radiation from the central AGN 



novae (e.g., 


ECohno et allboOl 


2007 


Kohno"2005'; 


Imanishi & Wada''2004l: Imanishi et al. 


2006a"b'; 


Imanishi & Nakanishi 


20061: 


Aalto et al.ii2007i: iMeiierink et al.ii2007i: iGarci'a-Burillo et al.H2007i: 


Loenen et al 


. boost lOarcfa-Burillo et al. 


20081: Perez-Beaupuits et all 20091). The strons UV and 



X-ray radiation from th e AGN and accretion disk co ul d affect both the dynamics and excitation of 
the rn olecular gas (e.g., Ohsuga & Umemura 200 1 al ibi: [Meijerink et al. 2007; van der Werf & et al 



2010|) . However, the radiation field from the AGN itself was not taken into account in the earlier 
estimates of molecular line emissions from hydrodynamical simulations. 

A prelimin ary estim ate of the potential effects of hard X-rays (E > \ ke V) on the molecular 
gas wa s done by IWPS09l using the X-ray Dissociated Region (XDR) models of lMeijerink & Spaans 
(120051). It was found that XDR chemistry may change the distribution of H2 around an AGN, if 
X-ray effects are explicitly included in the hydrodynamic model. The X-ray chemistry depends 
mainly on Hx/n, where Hx is the X-ray energy deposition rate and n is the number density of 
the gas (|Maloney et al.ll 19961) . Although the H2 abundance is robust in a clumpy medium like the 
one found in the hydrodynamical models, the temperature of the gas affected by an X-ray flux is 
expected to be a fa ctor of ~ 5 higher than that found in standard mo dels of a Photon Dorn inated 
Region (PDR; e.g. jHollenbach & Tielensll 19991) . for log(Hx)/n > 26 dMeijerink et al.ll2007l) . This 
is because the ionization heating by X-rays is more efficient than photo-electric emission by dust 
grains. 

Other molecular and atomic lines have also been suggested as tracers of t he AGN and starburst 
activit y in nearby galaxies {z < 1) as well as in high (z > 1) redshift galaxies. ISpaans & Meijerink 
(|2008|) studied the possibility of using '^CO and H2 emission lines to trace a young population 
of accreting massive (> 10^ Mq) black holes at redshifts z = 5 - 20 and radiating close to the 
Eddington limit. An enhancement in the intensities of various ^^CO transitions up the rotational 
ladder, as well as other molecular and atomic lines like '^CO, HCN, HCO^, [C i], [C 11], [O i] and 
[N ii], is also expected to be observed when X-ra y irradiation dominates the local gas chemistry 



(IMeijerink et al. 



20071 : ISpaans & Meijerinkll2008|) . Simulations of quasars at z ~ 6 with massive 



(10'^ - 10'^ Mq) halos and different merging histories showed that mid-7 '^CO lines are highly 
excited by a starburst, while high velocity peaks are exp ected to be produced by A GN-driven winds 
dNarayanan & et al.ll2008alJbl) . It was further found by Narayanan et al.l (|2009() that the compact 
^^CO spatial extents, broad linewidths and high excitation conditions observed in Submillimetre 
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Galaxies (SMGs) at z ~ 2 can be explained if SMGs are a transition phase of major merging events. 

In this work we use the XDR/PDR chemical model byl Meijerink & Spaans ( 2005 ) to estimate 
the abundances of more than 100 species (atoms and molecules) at each grid p oint in th e compu- 
tational box of the extended 3-D hydrodynamical models of an AGN torus by IWPS09l . We also 
estimate the actual X-ray flux emerging from the AGN, derived from the central black hole mass. 
Flux attenuation by photo-absorption of X-rays along the ray path and the distance from the central 
black hole is included. Thus, we estimate non-homogeneous abundances at each grid point that 
depend on the local density and impinging X-ray flux. A n extended version of the non-LTE 3-D ra- 
diative transfer code /33D byl Poelman & Spaans (20051) is used to compute the level populations o f 
any molecule or atom for which collision data exist in the LAMDaQ database (ISchoier et al. I l2005l) . 
Molecular and atomic line intensities and profiles are calculated with a line tracing approach for 
an arbitrary viewing angle. Model predictions for future ALMA observations of CO lines and the 
[C n] 158//m fine structure line are presented. The organization of this article is as follows. In 
Sec. |2]we describe the numerical method. The results and analysis are presented in Sec. |3l The 
final remarks are presented in Sec.|4l 



2. Numerical Method 



The three-dimensional hydrodynamic model of the AGN torus used in this work includes 
inhomogeneous density fields and mechanical heating effects due to turbulence and supernova 
explosions, with a resolution (pixel size) of 0.25 pc in diamete r. Detailed descriptions of the 
hydrodynamic equations and simulations can be found in lWPS09l The 3-D hydrodynamic model 
considers the formation and destruction of H2 in a self-consistent way, including formation of H2 
on grains, and the destruction of it by FUV radiation. This allows the code to compute the total 
density, local temperature (and velocity field) as well as the fraction of H2 at each grid element. The 
hydrody namical code runs for an equivalent time of > 3.5 Myr until it reaches a quasi-steady state 
(|WN02|) . and the gas forms a highly inhomogeneous and clumpy torus with some sp iral struc tures. 
It comprises a flared disk of H2 gas ~ 50 pc in diameter, and about 10 pc in height (|WPS09l their 
Fig.2c). 

One of the main criticism that hydrodynamical models receive in general, is that the density 
(and temperature) distribution that they show does not mimic c losely ac tual observations of the gas 
structure and distribution in galaxies and galaxy nuclei (like in lWPS09l their Fig.2a&b). However, 
observational data do not show the actual density or temperature of the gas either. The information 



' http://w w w. strw.leidenuniv.nl/ ~moldata/ 
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Fig. 1. — Number of grid elements in the 3-D hydrodynamical model that has both density n(H2) higher than 
10^ cm"-' and temperature lower than 10"* K, along the line-of-sight of the X-Y plane (left) and the Z-Y plane 
(right). These are the grid elements that would contribute the most to the emission of molecular gas irradiated by the 
X-ray flux emitted from the central SMBH. 



we get from observations is the intensity and distribution of particular atomic and molecular emis- 
sion (or absorption) lines, from which the ambient conditions (density, temperature and radiation 
field) can be estimated. Therefore, hydrodynamical models need to be complemented with atomic 
and gas chemistry that allow us to infer how the emission of different species would look like given 
the density and temperature structure obtained from the hydrodynamic simulations. 

One aspect that needs to be kept in mind is that not all th e grid elements shown in density 



and temperature distribution maps of the 3-D code by IWPS09I will contribute to the emission of, 
particularly, molecular lines. Figure [T] shows the number of grid elements with relatively cold 
temperature (T^ < 10^ K) and moderate density (n(H2) > 10^ cm~^) that can contribute to the 
molecular emission emerging along the line-of-sight of the face-on (X-Y plane) and edge-on (Z-Y 
plane) viewing angles. At higher temperatures and lower densities, the fractional abundances of 
molecular species would be very low (< 10~'") and their contribution to the molecular emission 
lines would be negligible due to the low coUisional excitation. The structure observed then is quite 
different than when considering the full density and temperature distribution at any cross section 
of the computational box. Although the maps of the main contributing grid elements represent 
a close estimate of the structure that we would expect to observe in molecular emission, they do 
not take into account the optical depth effects nor the (sub-)thermal excitation of the molecular 
energy levels that are treated in the radiative transfer calculations described in Sec. 12.31 The actual 
structure of emission lines then depends on the local density, temperature, and the radiation flux 
impinging at each grid element. 



The heating of the gas and dust by X-rays emanating from the AGN, as well as the chemical 
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abundances in the cold (< 500 K) gas, are not computed in hydrodynamical models. Including the 
(time dependent) chemical evolution at each step in the hydrodynamical simulations would take 
too long with the current computational resources. Theref ore, after a realization of th e state-of- 
art 3-D hydrodynamical simulation, the XDR/PDR code by lMeijerink & SpaansI (|2005|) is used to 
estimate the chemical abundances, based on the local density and impinging X-ray flux at each grid 
cell of the computational box. This code is depth dependent up to large columns (A^h ~ 10^^ cm~^), 
and considers a large (over 100 species) chemical network. In the following sections we describe 
the calculation of the X-ray flux, and the impact that it has on the chemistry and heating of the 
atomic and molecular gas around the AGN. 



2.1. The X-ray flux model 



For our 3-D hydrodynamical model we have a Mbh = 1.3 x 10^ M© SMBH (lWPS09h . so 
we can estimate the monochromatic luminosity of the AGN model at the rest frame wavelength 
/i = 510 nm as follows 



iL^(510nm) = 10^"^ x 



10"^ IM 



Mq 



l/h 



erg s" 



(1) 



With a = 5.7i:°;46 and b = 0.545 + 0.036 re.g.. iKaspi et al.ll2000h we have AL^iSlO n m) = 6.62 x 
lO^erg s"^ Using the same bolometric to monochromatic luminosity factor as in iKaspi et al 



(|200(^) we can determine the bolometric (total radiant energy) luminosity as L^o/ ^ 9AL^(510 nm) erg s ^ . 
The incident specific flux is assumed to have a spectral shape of the form 

FiiE) = Fo (y^) ^^P (-E/Ec) erg s'' cm'^ eY'\ (2) 

where E = hv eV, a ~ 0.9 is the characte r istic spectral index of the power-law compon ents of 
Seyfert 1 galaxies (e.g.. Pounds et all 1 19901 : iMadejski & et all 1 19951: Izdziarski et al.lll995l) . E^ is 
the f high energ y cut-off" which can be > 100 keV, 200 keV or 550 keV depending on the sample 
of AGNs (e.g.. IMadejski & et al.lll995 ). and Fq is a constant we estimate later to match the fraction 
of the total luminosity emitted in X-rays at the central grid point in the data cube. On the other 
hand, the lower energy cut-off" would depend on the shielding column density seen by the X-ray 
ffux at each grid point. Soft X-rays (photons with energy < 1 keV) cannot be eff"ectively attenuated 
by columns < 10^^ cm"^. However, as it is shown later in Sec.3.1, the column densities typically 
found in the inner ~ 25 pc of the torus are larger than 10^^ cm"^. Therefore, we only consider 
the hard X-rays between 1 keV and 100 keV as relevant for our X-ray chemical model. So we 
integrate eq.© over this energy range in order to obtain the hard X-ray ffux (Fi^rd) as 
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^100 keV , 
Jl keV ^ 



hard 



E \" 



1 keV; 



-E/E, 



dE 



erg s ' cm ^. 



(3) 



Considering that only ~10% (|Schleicher et alj 120101) of the total luminosity is emitted in X 



rays, we have that Fhard = 0.1 x Lhoi/Anr^, where ro is the distance from the central black hole 



which, for our purpose, is assumed to be the size of a grid cell (ro =|| ro(xo,yo, Zo) 
the central unresolved grid point. From this we find that Fq 1.4 x 10^ erg s~' cm 



1= 0.25 pc) for 
-2eV-i. 



For the rest of the cells, at position r{x, y, z) in the cube (a vector) the flux decreases not only 
with the square of the distance r =|| r(x,y,z) - ro{xo,yQ,Zo) II from the central black hole, but 
also because of the opacity r{E, r) of each grid cell at position r{x, y, z) along the radial path. The 
opacity is defined as 



T{E,f) = cTUE)Nn{r), 



(4) 



where apa{E) is the photoelectric absorption cross section per hydrogen nucleus, and N\i{r) is the 
total column density of hydrogen along the radial path from the central black hole to the position r 
in the computational box. The photoelectric absorption is calculated from all the species as 



(Tpa{E) = Y,^f'"'(Ti{E\ 



(5) 



with the total (gas and dust) elemental abundances, taken fromlMeiierink & SpaansI (|2005[) . 
and the X-ray absorption cross sections, cri(E), from Vemer & Yakovlev ( 1995I) . The total hard X- 
ray flux Fhardir) (erg s"' cm"^) impinging on an arbitrary grid cell at position r is then calculated 



as 



hard 



0.25 pc 



100 keV 



keV 



1 keV 



(6) 



Figure [2] shows the hard X-ray flux estimated for the 3-D hydrodynamical model by |WPS09 
in the X-Y plane 3.5 pc below the mid-plane of the AGN torus {left panel), as well as the flux in 
the actual X-Y mid-plane {middle panel). The right panel shows the flux in the Y-Z plane. The 
grid cells with diff'erent densities cause more or less absorption of the X-ray flux along the radial 
path from the central SMBH, producing shadow-like shapes and an inhomogeneous flux field. 

This total bolometric X-ray flux (from now on Fx) is used along with the total gas density 
of a grid cell as input parameters of the XDR/PDR chemistry code to estimate the abundances of 
several species. The formalism is described in the next section. 
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Fig. 2. — Impinging hard X-ray flux (in logarithmic scale and units of erg s"' cm"^) as seen in {left) the X-Y plane 
3.5 pc below the mid-plane of the AGN torus; {middle) the flux in the actual X-Y mid-plane; and {right) the flux seen 
in the Y-Z plane. Note that the X-ray flux distribution is not homogeneous. The shadow-like shapes are due to the 
X-ray absorption by the grid cells with different densities found along the radial path from the AGN. 

2.2. Chemical abundances and temperature 

For each grid point in the computational box we have the total gas density from the 3-D hy- 
drodynamic model. Since each grid point represents a physical (unresolved) scale of 0.25 pc, we 
also know the total column density that the impinging radiation flux will go through. Thus, the 
total gas density, the radiation flux, and the colu mn density of each grid po int are used as input 



parameters for the XDR/PDR chemical model by lMeijerink & SpaansI (|2005|) to compute the den 



sities and fractional abundances of different species at diff"erent column densities throughout the 
0.25 pc slab. In addition, we also get the temperature of the gas (as a function of the column den- 
sity) derived self-consistently from the chemical and thermal balance computed in the XDR/PDR 
code. Hence, we can compare the temperatures and H2 densities estimated from the X-ray free 
3-D hydrodynamical model with those computed considering the X-ray effects (see Sec. 13.41 ). 

Figure [3] shows the fractional abundances of some species as a function of the column density 
for impinging radiation fluxes of 160 erg cm~^ s"^ (left panel) and 1.6 erg cm~^ s"^ (right panel). 
The slab represents a single unresolved grid point in the computational box, with a fixed scale 
of = 0.25 pc and total hydrogen density hh = 10^ cm"^. This gives a total column density of 
A^H = n^x d ~ 10^^ cm^^, which is marked with a vertical slashed-dotted line. The XDR/PDR 
model, though, was executed to compute the abundances of the species up to a column of 10^^ cm~^ 
in order to show how deep the X-rays can penetrate depending on the strength of the radiation field. 
For the strong X-ray flux (left panel) and the actual column density (~ 10^^ cm"^) of the slab, the 
abundances of molecules like HCO^ will be negligible. Since the scale of the grid points is fixed, a 
slab denser than 10^ cm^^ would be required in order to observe a significant abundance of HCO^ 
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Fig. 3. — Left panel - Fractional abundances of various atomic and molecular species in a slab that represents one 
unresolved grid cell of d - 0.25 pc and density of ne = 10^ cm"^, in the 3-D hydrodynamic model. The X-ray flux of 
160 erg cm"^ s"' penetrates on the left side of the slab and affects the chemistry as it is absorbed throughout the column 
of gas and dust. The vertical dashed-dotted line indicates the total column density A^h = «h x ~ 8 x 10^^ cm"- of the 
slab. For this particular X-ray flux, a denser slab would be required to observe a significant abundance of molecules 
like HCO^. Right panel - Fractional abundance for the same slab as above, but with an impinging radiation flux of 
1.6 erg cm"^ s ' . Note the higher abundances of H2, CO, and HCO^at lower column densities, while the abundance 
of H, C and decrease earlier and faster throughout the slab. The overall temperature is also lower in this case. 

at larger (A^h > 10^^ cm"^) columns. Conversely, a weaker radiation flux of 1.6 erg cm"^ s"^ (right 
panel) impinging on the slab will produce a higher abundance of e.g., H2, CO, HCO"^ at lower 
column densities, while H, C and C^ would be less abundant, in a column averaged sense, than in 
the previous case. An overall lower temperature is also observed when the radiation flux is weaker, 
since X-ray photons are completely absorbed at a column of > 10^"^ cm"^. 

Since a grid point of the 3-D hydrodynamic model is unresolved, and the chemical abun- 
dances given by the XDR/PDR code depend on the column density, we estimate an abundance 
that is representative for the particular slab and for each species in the chemical network from the 
abundances observed throughout the slab of 0.25 pc. We compute the total fractional abundance 
< j?x > of the species Z='^CO, HCN, etc., as 

1 nx{l) dl 

< j^x >= (7) 

J Rh dl 

where nx(l) (cm"^) is the density of the species X at the layer I (cm) in the cloud, and is the 
total density of the slab, so the denominator is actually the total column density A^h (cm~^) of a 
particular grid point. 
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Because the variation of the abundance of a species through the slab is different for every 
species (e.g., [C i] is more abundant at the edge of the slab where T is high, while ^^CO is more 
abundant deep into the slab, where T is low), we require a gas temperature that is representative of 
the layers in the slab where the abundance of the species is higher (since those layers contribute the 
most to the line emission). Therefore, we compute an abundance- weighted average temperature, 
throughout an unresolved grid point, as 

/ ^x(l)T{l) dl 

< Tx >= , 8) 

/ ^xiD dl 

which gives different temperatures for different species in the same grid point. For instance, with 
the lower (1.6 erg cm"^ s"') X-ray flux {right panel in Fig. [3]) the total fractional abundance of CO 
is ~ 10""* with an abundance- weighted average temperature of ~ 37 K. While for C we have an 
abundance of ~ 3 x 10"^ and a representative temperature of ~ 66 K. 

Due to the large number of grid elements in the hydrodynamical model (256 x 256 x 128 ~ 
8.4 X 10^ data points) we require an optimized process in order to run the XDR/PDR code for 
the essential grid elements. Our criterion was to process only the grid points with total density 
larger than 100 cm"^ and temperature lower than 10"^ K. This is because lower gas densities will 
have little weight in the excitation of the molecular species, while higher temperatures would lead 
to mostly ionized gas. Once the abundances and abundance-weighted average temperatures have 
been determined for the selected grid elements of the 3-D hydrodynamical model, we can proceed 
to perform the radiative transfer calculations for the entire cube. 



2.3. 3-D radiative transfer and line tracing 



The advanced 3-D radiative transfer code ;S3D (|Poelman & Spaansll2005l) has been optimized 
for heavy memory usage, to be able to use the 256 x 25 6 x 128 elements data cube of the low 
resolution (0.25 pc scale) 3-D hydrodynamical models by IWPS09l . In principle, the temperature, 
density and velocity field derived from the hydrodynamical simulations can be used as the ambient 
conditions for the radiative transfer formalism. However, the temperature derived from the XDR 
model is found to be significantly different than the temperature obtained in the hydrodynamical 
model, as discussed in Sec. 13. 4[ A multi-zone radiative transfer approach is used, in which the 
calculation of the level populations in a grid cell depends on the level populat ions of all the other 



cells t hroug h the different escape probabilities connecting adjacent grid points (|Poelman & Spaans 

20051 bood) . 



The coUisional rates available in the LAMDA database (ISchoier et al.l 120051) are used in a 
similar way as in the 1-D and 2-D radiative transfer codes RADEX (|van der Tak et al.ll2007[) and 
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RATRAN (|Hogerheijde & van der Takll2000h . to calculate the level populations of different atomic 
and molecular species (e.g. [C i], [C ii], [O i], ^^CO, ^^CO, HCN, HCO+, HNC, CN, etc.). For 
this we use the commonly adopted main collision partner H2 for the radiative transfer calculations 
of all the molecules. Although the contribution of helium atoms to the total collision density for 
CO is just about 10"^, we also include (for completeness) He as an additional collision partn er by 
extrapolating (see AppendixjAj) the rate coefficients reported in lCecchi-Pestellini et al.l (120021) . For 
the case of [C 11] discussed in Sec. 13. 3[ we also use H and electrons as collision partners. With the 
exception of electrons, the density of all the collision p artners n (H) and n(H2) at each grid point 
is given by the hydrodynamical model (as described in (|WPS09|) ). since our aim is to know how 
the line emissions would look like for this particular model of an AGN torus. The hydrodynamical 
model, however, does not yet trace the evolution of electron density. So we use n(e~) derived from 
the XDR model. 

The line intensities, including kinematic structures in the gas, and optical depth effects, are 
computed with a ray-tracing approach for arbitrary rotation (viewing) angles about any of the three 
axes of the computational box. The emerging specific inte nsity is computed using the escape 
probability formalism described in lPoelman & SpaansI (|2005|) 



1 



(p{Vij)dz , 



(9) 



where dll has units of [erg cm s sr Hz^ ], n, is the population density in the i level, Aij the 
Einstein A coefficient, hvij the energy difference between the levels i and 7, S ij the source function 
for the corresponding transition, ll'^ivij) the local background radiation at the frequency of the 
transition, and t,^ is the cumulative optical depth that increases away from the observer, from the 
edge of the ray path to the z''^ layer. 



3. Analysis and results 

Once the level populations of particular transitions have been estimated with the radiative 
transfer code, and the line tracing at a particular inclination angle has been completed, the resulting 
two-dimensional emission map can be exported into a regular FITS data cube. 
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Fig. 4. — Top panels - Face-on view of the total column density (units of cm"^) A^h {left), the column density of 
molecular hydrogen A^(H2) (middle) and CO column N{CO) (right) in logarithmic scale. Bottom panels - Surface 
brightness maps of the / =1— >0toy = 9— >8 transitions of CO (in logio scale and units of erg s"' cm"^ sr"'), as 
observed at the surface of the face-on data cube. The emission of the lower CO 7=1— >0, 7 = 2-^1, and 7 = 3-^2 
transitions do not trace (or just with relatively fainter emission) the inner region of the torus, while the higher CO lines 
(from 7 = 4 — > 3 up to 7 = 9 — > 8) do trace the inner spiral structures, including the inner NLR. 



-14- 



3.1. CO maps 

The total hydrogen column density A^h, for a face-on viewing angle of the 3-D hydrodynami- 
cal model, is shown in the top left panel of Fig. |4l As expected, the total column density of the CO 
molecule N(CO) (top right panel) follows a similar distribution, although with about four orders of 
magnitude lower columns. The left and right bottom panels of Fig. |4] show the surface brightness of 
the CO 7=1^0 and 7 = 6 — > 5 lines, respectively. These correspond to the brightness observed 
at the surface of the face-on data cube (i.e., not scaled for an arbitrary distance to the source). 
Given the relatively lower upper energy state (~ 5.53 K) and critical density (~ 2 x 10^ cm"^ at 
100 K) of the CO 7=1^0 line, most of the warm and dense gas and structure is traced by 
the CO 7 = 6 — > 5 line instead. The emission of the lower CO 7 = 1 ^0, 7 = 2— > 1, and 
7 = 3 — > 2 transitions trace with relatively fainter emission (or not trace at all) the inner region 
of the torus, while the higher CO lines (from 7 = 4^3 up to 7 = 9^ 8) strongly trace the 
inner structures, including the inner Narrow Line Region (NLR) of the torus. This is also shown 
by the CO 3-2/1 - and 6-5/1 - line intensity ratios of Fig. [5] (top panels). This can be 
explained by an optically thick 7=1^0 line {N{CO) > 10'^ cm"^) and by the modest presence 
of cold (< 100 K) gas, specially at the inner +5 pc of the AGN torus. This indicates that, by just 
considering the excitation of CO, the 7=1^0 line will not always be a good tracer of hydrogen 
column density n^ in the central region (< 60 pc) of an AGN. The same can be concluded from the 
emission maps obtained considering an inclination angle of 45° about the X-axis. This is shown in 
Fig. Oof Appendix m 

As an exercise for comparison with future observations, we simulate a raster map of the AGN 
torus by adopting a distance D = 3.82 Mpc (the distance to NGC 4945) to the source, and by 
convolving the surface brightness maps with a single dish beam of FWHM=0.15" (about 1 1 pixels 
in the original map). This corresponds to a spatial scale of ~ 2.8 pc at the distance chosen, and 
gives a flux with units of 10"'^erg s"^ cm~^ after multiplying the surface brightness by the solid 
angle dO. = dR^/D^ subtended by the original pixel scale (dR = 0.25 pc) at the adopted distance of 
the source. A step size of one third the FWHM (~ 0.92 pc, or about 4 pixels) degrades the original 
image from 256x256 to a 61x61 pixels image. Figure [5] (bottom panels) shows the resulting flux 
maps of the CO 7 = 1 — > (left) and CO 7 = 6 — > 5 (right) lines. All the set of transitions 
from 7 = 1 — >0to7 = 9— >8is shown in Fig. [13] of the Appendix |B] The smearing effect of 
the relatively large beam produces the loss of the intricate structure observed in the original maps 
with 0.25 pc resolution shown in Fig. HI While a torus like shape becomes more evident in the CO 
7=1^0 line. 
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Fig. 5. — Top panels - Maps of the CO 3 - 2/1 - {left) and CO 6 - 5/1 - {right) line intensity ratios. Bottom 
panels - Face-on view of the flux (10~^^erg s"^ cm"^) of the CO 7=1^0 {left) and CO 7 = 6—^5 {right) lines, as 
mapped with a single dish beam of FWHM=0.15" 2.8 pc) and adopting a distance D = 3.82 Mpc to the source. 
Note how the relatively large beam smears out the intricate structure observed in the maps with the original resolution 
(0.25 pc) shown in Fig.|4l These lower resolution maps make more evident the fact that the 7=1^0 line traces the 
more diff'use outer gas while the / = 6 — > 5 traces the denser gas of the inner NLR. 
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3.2. The relation between CO and gas mass 



In the last few decades the integrated line fluxes of the lower ^^CO rotational lines (1-0, 2- 
1, a nd 3-2) have been u s ed to estimate the gas masses of molecular c louds in the Milky Way 
(e.g. lDickman et al.ll 19861 : ISolomon et al.ll 19871 : ISolomon & Barrettill99lL and references therein). 
These estimates hold for the Milky Way and nearby normal galaxies where the CO emission 
emerges from moderately dense (volume-averaged densities of n(H2) ~ 500 cm~^) giant molec- 
ular clouds (GMCs) in virial equilibrium (i.e., self-gravitating). For spherical clouds supported 
by isotropic random motions (e.g., turbulence) in virial equilibrium, the resulting inferred theo- 



retical conversion factors are X ~ 2 x 10 ^ ° [cm ^ (Kkm s ^) \ ] and a = 4.3 



(IStrong & Mattoxll 19961 : Dame et al.ll200ll : bickman et al.ll 19861 : ISolomon et al. 



Ma(Kkms-i)-i] 



1987h 



The virial approach of the optically thick CO lines can be extended to other gala xies, consid- 
ering an ensemble of virialized clouds, instead of a single one (IDickman et al.lll986h . In galactic 
nuclei and starburst galaxies, however, the assumption of an ensemble of individual gas clouds in 
virial equilibrium does not hold. In these environments the gas motions are due to a combination 
of gas and stellar mass co mponents, and the g a s is expected to be in a s moother configuration 
along a disk. Nevertheless J Pownes et al. ( 1993 ). Solomon et al. ( 1997 ). and Downes & Solomon 
(|1998h have shown that a slightly modified version of the CO luminosity to gas mass conversion 
factors can be derived in these environments. For luminous or ultraluminous infrared galaxies the 
inferred theoretical conversion factors range bet ween a = 0.8 and 1.6, or X = 3.7 - 7.3 x 10'^ 
dSolomon et al.ll 19971 : loownes & Solomonll 19981) . 



Since we know exactly what is the gas mass in our models, we can explore the behaviour of 
the CO-to-H2 conversion factor in our model of an AGN torus from the computed luminosities of 
several CO rotational lines. First we check the relation between the average density < > and 
the average abundance-weighted temperature of CO < Tqo > (computed through the line of sight 
or column) for each pixel of the map at the original resolution. We apply a similar criteria as in 
Sec. 2. 2 when computing the average density and temperature. That is, we use only the grid points 
with total densities larger than 100 cm~^ and temperatures Tqq lower than 5000 K to compute 
< Hh > and < Tqo >. This criteria, however, produces several pixels with < >= 0, specially 
in the outher region of the maps. Then we generate lower resolution raster maps by convolving 
the original resolution (0.25 pc) maps of < > and < Tqo > with difi'erent beam sizes (FWHM) 
corresponding to 1.75 pc, 4.75 pc, and 9.25 pc (with step or pixel sizes of ~FHWM/3), as it was 
done in Sec. 13.11 Those pixels with < hh >= are masked out in both density and temperature 
maps during the convolution process. From all the pixels of the raster maps we create scatter plots 
of < hh > versus < Tqo > at the difi'erent resolutions. 



The density and temperature maps, as well as the corresponding scatter plots, are shown in 
Fig.[6l The < hh > and < Tqo > maps (first and second column in Fig.[6l respectively) show denser 



<Tai> [K] 




<Tai> [K] 



Fig. 6. — Two dimensional maps (at the original resolution of 0.25 pc) of the average density < «h > (left panels), the 
average abundance-weighted temperature < Tco > of CO (middle panels), and scatter plot of < ne > versus < Tco >■ 
From top to bottom: same as in the top panels but convolving the maps with beam sizes (FWHM) equivalent to 1.75 
pc, 4.75 pc, and 9.25 pc, for a distance of 3.82 Mpc to the source. For the lowest resolution maps the scatter plot is 
shown with larger symbols for better visualization. All the scatter plots show the expected inverse relation between 
< hh > and < Tqo >'■ lower temperatures correspond to higher densities, and they converge to a linear relation as the 
beam size increases. 
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and colder gas in the center of the AGN torus, indicating an inverse relation between the average 
density and temperature (i.e., higher densities correspond to lower temperatures). The third column 
in Fig.[6]shows the scatter plots obtained considering all the pixels in the maps, where the expected 
inverse relation between < > and < Tco > is reproduced at all the resolutions. The relation 
between density and temperature becomes tighter and almost linear in the maps convolved with 
larger beams. 

Similarly, for each pixel of the CO maps we compute the total gas mass Mgas by adding the in- 
dividual masses of each grid point along the line of sight, and the CO luminosity L^^ derived from 
the CO intensity maps obtained in Sec. 13. II In this case we do not apply the criteria hh > 100 cm"^ 
and Tco < 5000 K as before, in order to obtain the actual total Mgas of the hydrodynamical model 
along the line of sight. However, we use the previously obtained average density map in order to 
mask the pixels with < >= in the gas mass and CO luminosity maps as well. This is because 
those pixels correspond to CO emission that is more than five orders of magnitude lower than the 
peak CO emission produced from our model. Given that is unlikely to have such large dynamical 
range in a real detector, the lower emission would be even below the noise level that could be 
found in maps obtained from real observations. Thus, not masking these pixels would cause an 
artificial bias toward lower luminosities during the convolution process. An alternative approach 
to masking pixels is to add random background noise to our intensity maps. However, the noise 
level to be added is arbitrary, so the results can be made at least comparable, and we are actually 
more interested in to show what is really coming out from the hydrodynamical and XDR models. 



-1 



Figure |7] shows the scatter plots of the gas masses (MgaJ and the CO luminosity (L^q [Kkm s 
pc^]) of the 7=1— >0, 7 = 6— >5 and 7 = 9 — > 8 transitions as obtained from the low 
resolution (0.25 pc) hydrodynamical model. We consider a 36% mass correction to account for 
the mass of helium atoms. The straight line corresponds to the gas mas s estimated assuming 



a luminosity-to-gas mass conversion factor a of Mgi,s/L'^o = or = 0.8 (e.g. ISolomon et al.lll997 



Downes & Solomonlll998 ). The result of the original resolution shows a large scatter of gas 



masses 

for a given CO luminosity. When considering the source at a nominal distance of 3.82 Mpc, 
however, the beam averaged masses and luminosities result in a tighter correlation, but still with a 
considerable scatter. The relation between L^q and Mgas is clustered just at the higher luminosity 
and mass ranges as the beam size increases. That is, the pixels with lower mass and luminosity 
are missing (beam smeared) in the lower resolution maps. The higher the CO 7-line and the 
larger the beam used, the closer is the relation found with our models to the luminosity-to-gas 
mass conversion factor proposed in the literature. The a factor was estimated based on the lower 
transitions (7= 1 — >0, 7 = 2— > 1,7 = 3^2), which are way off" in our models because 
there is very little cold (Tco < 100 K) gas in the inner 60 pc. Besides, the relatively large scatter 
observed at all resolutions indicates that many of the clumps where the CO emission emerges from 
are not in virial equilibrium. On the other hand, the CO 7 = 6 — > 5 and 7 = 9^8 transitions 
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Fig. 7. — Scatter plots between the ensemble of gas masses (Mgas [Mq]) and the CO luminosity (L^q [Kkm s"' pc^]) 
of the (from left to right) 7=1— >0, 7 = 6— >5, and 7 = 9 — > 8 transitions derived from the low resolution (0.25 pc) 
hydrodynamical model, and convolved (beam averaged) with beam sizes (FWHM) equivalent to 0.25 pc (the original 
scale of the model), 1.75 pc, 4.75 pc, and 9.25 pc (from top to bottom, respectively), for a distance of 3.82 Mpc to the 
source. The straight line corresponds to the gas mass estimated assuming a luminosity-to-gas mass conversion factor 
a of Mgas/L[~,Q - a - 0.8 (e.g., Solomon et al. 1997; Downes & Solomon 1998). The gas mass used here includes a 
36% correction to account for helium. Note the different scales of the axis for different resolutions and 7-lines. 
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show similar trends. This implies that the higher-7 line does not provide significant additional 
information compared to the 7 = 6 — > 5 line. 

We can confirm then a tight correlation between the beam averaged luminosity of the higher 
CO 7-lines (from 7 = 5— >4 to 7 = 9^ 8) and the gas mass of the AGN torus. However, the 
larger scatter seen when using higher resolution (smaller beam sizes) should be taken as a warning 
sign for future higher resolution observations like ALMA, which may no longer show a linear 
correlation between CO luminosity and gas mass, or may introduce a larger range of correlations 
for ensembles of clouds with different density and high temperature structures. So, using only 
mid-7 CO lines (like 7 = 6 ^ 5) for gas mass determinations in AGN tori would be the most 
reliable approach. 



3.3. The [C ii] ISS/um fine structure line 



A considerable amount of gas close to the SMBH is predicted to be at very high temperatures 
(T/, > 1000 K) in the 3-D hydrodynamical model. At these temperatures hydrogen is found mostly 
in atomic form, and other atoms like carbon are mostly ionized. Therefore, in this section we 
present the results of the radiative transfer calculations done for [C ii]. We use the collision rate 
coefficients of [C ii] from the LAMDA database. In contrast with the CO molecule, we considered 
not only the molecular hydrogen as colli sion partner, but also the atomic hydrogen and the elec- 
trons, using the collis ion data reported by iFlower & Launayl (|l977|) . lLaunay & Rouefj (|l977[) . and 
Wilson & Bell tooi) . These are the most relevant collision partners for [C ii] in an AGN envi- 
ronment, since their densities are expected to be higher than ^(Hi) in the very hot gas (Tt > 1000 
K) regions. As described in Sees. 12.21 and |2.3[ the fractional abundance of the electrons is derived 
from the XDR mode l. Whilst the den sity n(H) is computed from the hydrodynamical model as 
n(H) = nu- 2n(H2) (iWada et aklboogl) . 



Figure [8] shows the face-on view {left panels) of the total column density A'^(e~) (cm^^) of 
electrons, the surface brightness (in logiQ scale and units of erg s"' cm"^ sr"') of the [C ii] 158 /im 
emission at the original spatial resolution (0.25 pc), and at ~ 2.8 pc resolution (for a distance D = 
3.82 Mpc to the source) after convoling the original map with a single dish beam of FWHM=0.15". 
This results in a flux with units of 10"'^erg s ' cm"^ after multiplying the surface brightness by the 
solid angle as described in Sec. 13.11 Comparing with the CO emission shown in Figs. |4] and |5l it 
is clear that the [C ii] 158 yum emission traces mostly the central region (inner NLR) of the AGN 
torus, and it is a better tracer of the hot regions than the mid-7 CO lines. 

The right panels of Fig. [8] shows the same as in the left panels, but with an inclination angle 
of 45" about the X-axis. The viewing angle produces a slightly larger column of [C ii] observed 
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through the line of sight, which in turn results in a ~ 30% brighter [C n] 158 jum peak emission. 
The [C ii] 158 /im emission (rest frequency ~ 1900 GHz) from galaxies at redshift z > 1 will be 
observable with ALMA. Unfortunately, we cannot obtain high resolution maps of [C ii] with the 
current specifications of ALMA. Even if we consider both the largest baseline of 16 km and the 
highest frequency of 950 GHz (band 10) that ALMA is expected to have in the future, we could 
get a resolution of FWH]V@~60/16(km)/950(GHz)~0.004". This in turn gives a spatial scale of 
32 pc for a source at z=l (with an equivalent angular distance of 1658.6 Mpc, in a LambdaCDM 
cosmology with Ho=71 km/s/Mpc). That is, the whole region we show in Fig.[8]would correspond 
to just one or two pixels in the ALMA maps of sources at z = 1. 



3.4. Temperature and density driven by X-rays 

In eql6]we assumed that the X-ray flux is spherically symmetric with respect to the central 
SMBH. If we assume instead that the X-rays are emitted in a preferential direction perpendicular 
to the mid-XY-plane of the accretion disk (and of the AGN torus as a whole), we can consider the 
X-ray flux emerging from a Lambertian object. That is, the radiation flux impinging on each grid 
point of the cube is proportional to the cosine of the viewing angle 6, with respect to the vertical 
Z-axis. This means that the X-ray flux would be negligible in the disk (mid-plane) of the AGN 
torus. 

In order to compare the impact that the two diff'erent X-ray flux distributions have on the 
temperature and molecular hydrogen density of the gas, we compare the temperature Thyd obtained 
from the hydrodynamical model with the H2 abundance-weighted average temperature Txdr (eqU]) 
derived from the XDR chemical model using both the spherical and Lambertian X-ray fluxes. We 
take a strip volume of 64 x 1.25 x 1.5 pc^ along the X-axis, and around the center of the Y- 
axis (AY = 0) and Z-axis (AZ = 0) of the 3-D cube. We use this thin volume so we can have 
similar X-ray fluxes (decreasing mostly with radial distance) impinging at each grid element of 
the 1.25 X 1.5 pc^ slices of the volume. We computed the average temperature and H2 density of 
the 1.25 X 1.5 pc^ slices at each AX grid element. At the resolution of 0.25 pc/element we have 
30 grid elements per slice, which is a good compromise between a representative number of grid 
elements, and a fairly constant impinging X-ray flux at each slice. 

The top panel of Fig. [9] shows the average temperature Thyd (K) estimated in the 3-D hy- 
drodynamical model {solid line) at each AX grid element. The corresponding average (of the H2 
abundance-weighted average) temperature Txdr (K) obtained from the XDR chemical model is 



■^http://science. nrao.edu/alma/specifications.shtml 



-22- 




Fig. 8. — Left panels - Face-on view (from top to bottom) of the total column density N{e^) (cm"^) of electrons in 
logarithmic scale, the face-on view of the surface brightness (in logio scale and units of erg s"' cm"^ sr"') of the [C ll] 
158 fim emission, and the same map as above but convolved with a single dish beam of FWHM=0.15" (~ 2.8 pc at an 
adopted distance D - 3.82 Mpc to the source), which gives a flux in units of 10"''' erg s"' cm"^. Right panels - Same 
as in the left panels, but with an inclination (viewing) angle of 45° about the X-axis. 
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Fig. 9. — Top left - Average temperature Tuyd (in units of K) of the gas along a 64 x 1.25 x 1.5 pc-' strip volume, 
as estimated in the 3-D hydrodynamical model {black line) and the corresponding H2 abundance-weighted average 
temperature Txdr obtained from the XDR chemical model (gray line) using the spherical X-ray flux (eq|6ll. The filled 
circles show the actual data points obtained with the XDR model in grid cells with Thyd < 10^ K. Top right - Average 
Txdr/Thyd ratio (gray line + filled circles) and the average spherical X-ray flux (solid black line) Fx (erg s"' cm"^) 
in logio scale. The standard deviation at each AX off'set is shown by the error bars. The relative average temperature 
is directly related to the impinging flux at each grid point. The average temperature Txdr is predominantly higher 
than Thyd in the inner 20 pc around the center of the AGN torus. The dashed line indicates where TxdrIThyd - 1- 
The Bottom panels show the same as above, but for the Lambertian X-ray flux. The average temperature Txdr is still 
higher than Thyd, but the difl'erence is smaller than for the spherical X-ray flux, particularly beyond +5 pc from the 
central SMBH. 
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Fig. 10. — Top left - Average density of molecular hydrogen n(H2)HYD (in units of cm"^) of the gas along a 64 x 
1 .25 X 1 .5 pc^ strip volume, as estimated in the 3-D hydrodynamical model {solid line) and the corresponding average 
n(H2)xDR density obtained from the XDR chemical model {gray line + filled circles) using the spherical X-ray flux 
(eq|6l). The filled circles show the actual data points obtained with the XDR model in grid cells with Thyd < 10"^ K. 
Top right - Average n(H2)xDR/n(H2)HYD ratio {gray line + filled circles) and the average spherical X-ray flux {solid 
black line) Fx (erg s"' cm"^) in log\Q scale. The standard deviation of the relative density at each AX off'set is shown 
by the error bars. The dashed line shows where n(H2)xDR/n(H2)HYD = 1 ■ The Bottom panels show the same as above, 
but for the Lambertian X-ray flux. The turnover in the relation of the hydrogen density is still at about AX = 10 pc 
from the central SMBH, but «(H2)xdr can be up to 10^ times higher than for the spherical X-ray flux beyond +10 pc. 
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shown with a gray line. Only the grid elements with Thyd < 10"^ K were used in the XDR model, 
and they are shown with filled circles. The average I'xdr/T'hyd ratio (gray line + filled circles) and 
the average X-ray flux (solid black line) Fx (erg s"' cm"^) are shown in the bottom panel of Fig. |9l 
The average relative temperature is directly related to the impinging flux at each AX grid element, 
and decreases as the X-ray flux decreases. The temperature derived from the XDR model is higher 
(T'xdr/T'hyd > 1) than the one estimated in the hydrodynamical model. Hence, the presence of 
X-rays has an undeniable eff"ect on the thermodynamics of the AGN torus, up to at least 60 pc. 

In Fig. [TO] the average density of molecular hydrogen n(H2)HYD along a 64 x 1.25 x 1.5 pc^ 
strip volume is shown. The density computed for the 3-D hydrodynamical model is shown by 
the solid line, and the corresponding average n(H2)xDR density obtained from the XDR chemical 
model, using the spherical (top panels) and the Lambertian (bottom panels) X-ray fluxes, is shown 
with the gray line + filled circles. Only the data points for grid cells with Thyd < lO'^ K and 
«(H2)hyd > 10"^ cm"^ are shown in the figure. 

The relative H2 density seems to be inversely related to the impinging flux. That is, the average 
n(H2)xDR density derived from the XDR model is lower (by factors up to ~ lO'^) than the average 
H2 density of the hydrodynamical model. This is observed in the inner ±10 pc region around the 
center of the AGN torus. This is mostly the consequence of relatively thin slabs (Ny{ < 10^^ cm~^) 
being irradiated by a rather strong (Fx > 1.6 erg s"^ cm"^) X-ray flux in the proximity of the torus 
center, as described in Sec. 12.21 and shown in FigjSl However, beyond AX ~ 10 pc the XDR H2 
abundance exceeds the one in the pure hydrodynamical model. This change is explained by the 
background star formation considered in the hydrodynamical model at those distances. The FUV 
is a more efficient destroyer of the H2 gas than X-rays (Meijerink & Spaans 2005). 

With a spherical X-ray flux n(H2)xDR can be up to lO'^ times higher than n(H2)HYD, while it 
can be about 10^ higher if a Lambertian X-ray flux is considered instead. In the inner region (AX ~ 
-2 pc.) a density «(H2)xdr a few times higher than n(H2)HYD is observed. This is the consequence 
of a weaker X-ray flux, with respect to the spherical radiation flu x. This is consistent with the 



fluctuations of the column density distribution of H2 explored by IWPS09I for diff'erent viewing 
angles and, as expected, the largest A'^(H2) columns are found at a viewing angle ~ degrees (i.e., 
edge-on). These facts imply that molecules will tend to disappear in the central (< 10 pc) region. 
But, depending on the viewing angle, and for total hydrogen columns > lO^'^ cm~^ like in the case 
of, for instance, the LIRG NGC 4945, molecules can survive and emission lines like, e.g., high-7 
CO, [C 11], [Ne 11] and [Ne v], can be bright. In all, we conclude that the H2 abundance in the 
AGN torus is strongly aff"ected by the black hole (< 10 pc) and star formation (> 10 pc) (see also 



Schleicher etal.ll2010 ). 



We note, though, that the nature of the models used here to derive temperature and density of 
the gas (a static X-ray driven chemical model and an hydrodynamical X-ray free model) are differ- 
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ent, and a comparison between their corresponding derived temperatures and densities is merely 
intended to motivate the need for a joint XDR-hydrodyna mical model for the therm odynamics of 
AGN tori. A first attempt to this effect has been made by iHocuk & SpaansI (|2010|) for individual 
~1 pc molecular clouds close to a SMBH. 



4. Final remarks 

We compared the total hydrogen column density, N(CO) and CO 7 = 1 —>0 to 7 = 9^8 
line intensities, and found that the mid-7 CO lines are excellent probes of density and dynamics, 
but the low-7 CO lines are not good tracers of in the central (< 60 pc) region of the AGN torus. 
The analysis of the Xco/o: conversion factors indicated that only the higher-7 CO lines will show a 
linear correlation with the gas mass in AGN tori at lower spatial resolutions (~ 9 pc). But at higher 
resolution (< 5 pc) different proportionality factors (or no correlation at all) pertained between the 
CO lines and the total gas mass in AGN tori. We also determined that the [C ii] 158 fim emission 
will trace mostly the central region of AGN tori, detectable (but not resolved) by ALMA in z > 1 
galaxy nuclei. 

We found that the presence of X-rays has an undeniable effect on the thermodynamics of the 
AGN torus, up to at least 60 pc. An important implication of this is that circum-nuclear star forma- 
tion could be suppressed in the central ~5 pc. This can shed light on the starburst- AGN connection. 



Self-c onsistent UV/X-ray radiation-chemical-hydrodynamical simulations (e.g. iHocuk & Spaans 



2010|) will allow us to explore this theoretically, and their predictions can be confirmed (and used 



for data interpretation) by ALMA in the near future. 

With a rest frequency of 691.5 GHz the CO 7 = 6^5 line can be observed with the ALMA 
band 9 receivers, which will be the highest frequency band available when the early science begins 
with > 16 antennas. Considering a minimum baseline of 250 m for the compact configuration, 
we would have an angular resolutior^ of FWHM«0.35" that will allow us to resolve structures 
of ~7 pc at a nearby distance of 4 Mpc (about the distance to NGC 4945) and of ~25 pc at a 
distance of 15 Mpc (roughly the distance to NGC 1068). In the near future, however, the higher 
sensitivity (with >50 antennas) and the availability of longer baselines of up to ~15 km will provide 
angular resolutions of FWHM«0.006" with ALMA band 9, which will allow the study of structures 
between ~0.1 pc and ~0.4 pc, respectively, at the distances mentioned above. Therefore, the spatial 
scales (>0.25 pc) that we probe with our simulations match the angular resolutions provided by 
ALMA. 



^http://science. nrao.edu/alma/earlyscience.shtml 
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A. Rotational excitation of CO by He 



We used the rate coefficients for pu r e rota tional de-excitation of CO by collisions with He 
atoms reported in ICecchi-Pestellini et al.l (l2002h . The original rate coefficients are given for the 
first 15 rotational levels and for ten different temperatures from 5 to 500 K. In order to extend the 
available rates to higher rotational levels and temperatures, we followed the methodology for linear 
molecules described by lSchoier et al.l (|2005|) . which was used to produce the LAMDA database. 



We first extrapolated the downward coUisional rate coefficients (A7 = 7„ — > 7,, 7„ > 7/ ) 
in temperature (up t o 2000 K) using t he modified version of the analytic approximation given by 
de Jong et al.1 (|l975|) and presented by lBieging et al.1 (Il998|) : 



r„/ = A(A7) y exp [-5(A7) y"^] x exp [-C(A7) y'"-] , 



(Al) 



where y = ^E uilkT and the three pararneters A, B, and C are determined by least-squares fits to the 
original set of ICecchi-Pestellini et al.l (|2002h rate coefficients for each A7. Then we extrapolated 
the coUisional rate coefficients to include higher rotational levels (up to 7 = 40) by fitting the rate 
coefficients (in natural logarithmic scale) connecting the ground rotational state to a second order 
polynomial 



l^iyjo) = a + bJ + cJ , 



(A2) 



with a, b, and c parameters deterr nined from the fit, fo r each temperature. The Infinite Order 
Sudden (lOS) approximation (e.g.. iGoldflam et al.lll977h was used to calculate the whole matrix 
of state-to-state rate coefficients from the coefficients connecting the ground state jlo 
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Fig. 1 1 . — Top left panel - Collisional rate coefficients (cm^ s ') for CO with para-H2 as collision partner Top right 
panel - Collisional rate coefficients for CO colliding with He with the wrong Wigner 3-j function. Bottom panel - 
CoUisional rate coefficients for CO-He collision partners with the correct Wigner 3-j function. 



L=J+J' 



L=\J-J'\ 



r..=(27' + l) 2 (2L+l)(^ ^ 



no, 



(A3) 



where the term 



J' 




L 




(A4) 



is the Wigner 3- j symbol that designates the Clebsch-Gordan coefficients (e.g. jTuzun et al.lll998|) . 
and references therein). The lOS approximation provides an accurate description of the collisional 
rates if the rotational energy differences are small compared to the kinetic energy of the colliding 
molecules. In cases where this condition is not satisfied, it is possible to approximately correct 
for the d eviations by multiply ing t he summation in eq.(IA3|) with the adiabaticity correction factor 
given bv bepristo et all (11972) and lMcKee et all (Il982h 



A{L,J)= ' with a = 0.13 5o/ ^ , (A5) 

6 + {ajy \T I 

where Bq is the rotational constant of the colliding molecule in cm^^ {Bq = 1.9225 cm"^ for CO), 
Z = 3 A is a typical scattering length, is the reduced mass of the colliding system in amu (/i « 3.5 
amu for CO-He), and T is the kinetic temperature. 

However, the A(L, J) correction factor should be used only if Ei> Ej and {Ei - Ej) ^ Ek, 
where E jl is the energy of the CO rotational levels L, J and E^ is the kinetic energy of the collision 
partners. The top left panel of Fig[TT] shows the deviations introduced by the A{L, J) factor when 
used arbitrarily to extrapolate the rate coefficients of CO colliding with He. The top right panel of 
FigHU shows similar discontinuities in the extrapolated rate coefficients between CO and para-H2 
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presented in the current LAMDA molecular data. Similar deviations are observed for ortho-H2. 
This means that the extrapolated LAMDA molecu lar data for CO ne ed to be corrected. A lthough, 
the original CO-H2 rate coefficients obtained from lFloweii (1200 1|) and lWemli et al.l (|2006|) go up to 
J = 29, and since we do not explore CO transitions above 7 = 20 we can still use the LAMDA 
molecular data without corrections. 

In the case of the CO-He colliding system the conditions for using the A(L, J) adiabaticity 
factor are not satisfied for the temperatures and energy levels consid ered here, and the IP S ap- 
proximation given by eq. (|A3l ) yields results with 10%-15% accuracy (IGoldflam et al.lll977|) . The 
final extrapolated rate coefficients used in this work for the system CO-He are shown in the bottom 
panel of FigfTTl 



B. CO maps at 45" inclination 
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Fig. 12. — Top panels - Maps with a 45° inclination about the X-axis of the total column density (units of cm"^) 
A^H (^e/f), the column density of molecular hydrogen A^(H2) (middle) and the CO column N{CO) (right) in logarithmic 
scale. Bottom panels - Surface brightness maps of the 7=1— »0to7 = 9— >8 transitions of CO (in log lo scale and 
units of erg s"' cm"^ sr"'), as observed at the surface of the face-on data cube. The larger columns seen through the 
line of sight with the 45" inclination produce higher emissions of the CO transitions with respect to the face-on maps 
of Fig. m However, the higher CO lines (from 7 = 4 — » 3 up to 7 = 9 — > 8) are stiU better tracers of the inner region 
of the torus. 
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Fig. 13. — Maps of the flux (erg s ' cm"^) of the CO transitions (from y=l->0toy = 9-»8) with a 45° 
inclination about the X-axis and convolved with a single dish beam of FWHM=0.15" (~ 2.8 pc). We adopt a distance 
D = 3.82 Mpc to the source. Although the J = 3 2 transition does not trace the warmer gas of the inner NLR, it is 
the brightest emission line of this configuration (see the colour scale). 
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